Identification of Rickettsia spp. in Ticks Removed from Tick-Bitten Humans in Northwestern Spain

Simple Summary We conducted a tick surveillance study in northwestern Spain with the aim of identifying tick-associated Rickettsia species among ticks removed from humans and to track how tick populations and their associated pathogens have changed over the years. Five genera and eight species of ticks were collected, with Ixodes ricinus being the most frequently found. By comparing our data with previous reports, a clear geographic and seasonal expansion of tick populations and the associated Rickettsia species is observed, indicating that the spatiotemporal patterns of ticks and bacteria have changed over the years. Considering the medical relevance of tick-borne rickettsioses, evaluating the infection risk to humans in tick-infested areas, as well as detecting spreading species, is essential for establishing adequate prevention and control measures. Abstract Tick-borne rickettsioses (TBRs) are distributed worldwide and are recognized as important emerging vector-borne zoonotic diseases in Europe. The aim of this study was to identify tick-associated Rickettsia among ticks removed from humans, and to track how tick populations and their associated pathogens have changed over the years. For this purpose, we conducted a tick surveillance study in northwestern Spain between 2018 and 2022. Ticks were morphologically identified and analyzed for the presence of rickettsial pathogens through the amplification of the citrate synthase (gltA) and the outer membrane protein A (ompA) genes. PCR products were sequenced and subjected to phylogenetic analyses. We collected 7397 ticks, with Ixodes ricinus being the species most frequently isolated. Based on the PCR results, Rickettsia DNA was detected in 1177 (15.91%) ticks, and 10 members of Rickettsia were identified: R. aeschlimannii, R. conorii subsp. conorii, R. conorii subsp. raoultii, R. massiliae, R. monacensis, R. sibirica subsp. mongolitimonae, R. slovaca, R. helvetica, Candidatus R. barbariae, and Candidatus R. rioja. Some of these Rickettsia have gone previously undetected in the study region. There is clear geographic and seasonal expansion not only of tick populations, but also of the associated Rickettsia. The comparison of our data with those obtained years ago provides a clear idea of how the spatiotemporal distributions of ticks and their associated Rickettsiae have changed over the years.


Introduction
Rickettsioses are worldwide zoonoses caused by intracellular Gram-negative bacteria of the genera Rickettsia and Orientia (order Rickettsiales, family Rickettsiaceae).These infectious diseases are transmitted by different arthropod vectors, such as lice, fleas, mites, and ticks, the latter being involved in the transmission of most pathogenic Rickettsiae to Insects 2024, 15, 571 2 of 19 humans.Tick-borne rickettsioses (TBRs) are distributed worldwide and are recognized as important emerging vector-borne zoonotic diseases in Europe [1,2], with cases being exclusively due to members of Rickettsia [3,4].
The genus Rickettsia currently comprises more than 30 species of intracellular coccobacilli bacteria, validly published under the International Code of Nomenclature of Prokaryotes (ICNP), and almost 90 candidate species (https://lpsn.dsmz.de/genus/Rickettsia;accessed on 20 April 2024).Among them, half are human pathogens [5], and eight are transmitted by different species of ticks in Europe [6].Both pathogenic and nonpathogenic Rickettsia species have been classified into five distinct well-supported clades: (i) two spotted fever groups (SFGI and SFGII, also known as transitional groups (TRGs)); (ii) the typhus group (TG); (iii) the basal bellii group (BG); (iv) the canadensis group (CG); and (v) the Tamurae/Ixodes Group (TIG) [7][8][9].SFGs are the most numerous group, with more than 56 validated and candidate species [10], although this number will likely change over time as novel SFG Rickettsia species continue to be identified.
In Europe, most SFG Rickettsiae are transmitted by ticks of the family Ixodidae.Specifically, the genera Dermacentor, Haemaphysalis, Hyalomma, Ixodes, and Rhipicephalus are considered the principal sources of infections naturally transmitted by ticks [11].According to the European Centre for Disease Prevention and Control (ECDC), most cases of rickettsiosis have been reported in Italy, Portugal, and Spain.
Focusing on Castilla y León (northwestern Spain), the second largest regions in Europe, with a wide variety of ecosystems and bioclimatic conditions, there is a recent study on human-biting tick species and their spatial and temporal distributions [13], but their associated rickettsial pathogens were not analyzed in this study.Indeed, as far as we are aware, there is only an exhaustive study on ticks removed from people in Castilla y León and their associated rickettsial pathogens by members of our group, but it was carried out in the period 1997-2002 [14].Since then, no additional studies have been conducted.The observed changes in tick populations by Vieira Lista et al. [13] make it essential to re-examine the current situation of pathogens transmitted by tick populations in this area.
The main purpose of this study was to fill this knowledge gap by identifying tickassociated Rickettsia species in ticks recently removed from humans in Castilla y León and comparing them with those found in the same area in the past.By doing so, we will be able to examine the spatiotemporal changes in the presence of ticks and their pathogens.Ultimately, this will allow us to evaluate the risk of TBR in tick-infested areas, as well as to detect spreading species and emerging diseases, which is essential for revealing its epidemiology and establishing adequate prevention and control measures.

Tick Collection and Identification
Ticks were collected over a 5-year period (2018-2022) from people who went to Primary Health Care Centres and Hospital Emergency Services in Castilla y León (northwestern Spain) (Figure 1), for their removal through a program for the prevention and control of tick-borne anthropozoonoses.
For each tick received, other data, such as sex and degree of feeding, as well epidemiological characteristics of the patients (age, sex, geographic location at the t the bite(s), occupation, and the anatomical location of the tick bite), were collected.
All ticks examined are kept at −20 °C in the tick collection of the Laboratory of tious and Tropical Diseases of the University of Salamanca (Salamanca, Spain).

DNA Extraction, Amplification, Purification and Sequencing
Genomic DNA was extracted from ticks using a NucleoSpin ® Blood Kit (Mac Nagel GmbH & Co. KG, Düren, Germany) according to the manufacturer's instru DNA was individually extracted from adult ticks.Only when several immature tick removed from the same patient, a DNA pool was made (specifically, six pools wer pared with 15 immature ticks removed from six people).DNA extractions were sto −20 °C until use.
To detect the presence of rickettsial genetic material in the DNA samples ob from the ticks analyzed, we used the citrate synthase gene (gltA) and the outer mem protein A gene (ompA), which are present in all members of the SFG but not in clos tives [19].Both genes were targeted by conventional PCR, including negative and po controls in each assay, using previously reported primers [20,21] and protocols (Tab and S2).Amplification reactions were performed in a Biometra ® T-personal Thermo (Whatman Biometra, Göttingen, Germany).
PCR products were visualized on a 1.5% agarose gel under a UVP Biodoc-It ® 2 ing system (Analytik Jena, Jena, Germany).All positive PCR products were purified a Nucleospin ® Gel and PCR Clean-up kit (Macherey-Nagel GmbH & Co. KG, Düren many) and quantified with a NanoDrop ND-1000 Spectrophotometer (Nanodrop nologies, Wilmington, NC, USA) with ND-1000 V3.8.1 software.Once purifie Ticks were sent to the Laboratory of Infectious and Tropical Diseases (e-INTRO), Faculty of Pharmacy at the University of Salamanca (Spain), for identification and pathogen detection through molecular analyses.At their arrival, each tick was morphologically identified under a binocular lens according to life stage and sex using taxonomic keys and reference works [15][16][17].Each specimen was identified to the species level, except for those ticks corresponding to Rh. sanguineus sensu lato, which were identified only to the group level due to the taxonomic issues of this species [18].
For each tick received, other data, such as sex and degree of feeding, as well as the epidemiological characteristics of the patients (age, sex, geographic location at the time of the bite(s), occupation, and the anatomical location of the tick bite), were collected.
All ticks examined are kept at −20 • C in the tick collection of the Laboratory of Infectious and Tropical Diseases of the University of Salamanca (Salamanca, Spain).

DNA Extraction, Amplification, Purification and Sequencing
Genomic DNA was extracted from ticks using a NucleoSpin ® Blood Kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) according to the manufacturer's instructions.DNA was individually extracted from adult ticks.Only when several immature ticks were removed from the same patient, a DNA pool was made (specifically, six pools were prepared with 15 immature ticks removed from six people).DNA extractions were stored at −20 • C until use.
To detect the presence of rickettsial genetic material in the DNA samples obtained from the ticks analyzed, we used the citrate synthase gene (gltA) and the outer membrane protein A gene (ompA), which are present in all members of the SFG but not in close relatives [19].Both genes were targeted by conventional PCR, including negative and positive controls in each assay, using previously reported primers [20,21] and protocols (Tables S1 and S2).Amplification reactions were performed in a Biometra ® T-personal Thermocycler (Whatman Biometra, Göttingen, Germany).
PCR products were visualized on a 1.5% agarose gel under a UVP Biodoc-It ® 2 imaging system (Analytik Jena, Jena, Germany).All positive PCR products were purified using a Nucleospin ® Gel and PCR Clean-up kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) and quantified with a NanoDrop ND-1000 Spectrophotometer (Nanodrop Technologies, Wilmington, NC, USA) with ND-1000 V3.8.1 software.Once purified, the amplicons were Sanger sequenced in one direction at the Sequencing Service of the University of Salamanca, Nucleus (Salamanca, Spain), with the same primer sets used for amplification.

Sequence Edition, Alignment and Phylogenetic Analyses
Sequence editing, which involved primer region trimming and consensus assembly, was performed in Geneious v7.1.9.(Biomatters, Ltd, Auckland, New Zealand; www.geneious.com).Edited sequences were subjected to BLAST searches against the NCBI GenBank database to verify that they did not correspond to ticks, and the sequences were assigned a putative identity based on the BLAST results (Table S3).Sequences that were too short relative to the global length of the alignment were not considered for subsequent phylogenetic analyses.All 581 new sequences used for phylogenetic purposes were submitted to GenBank under accession numbers PP552066-PP552646 (Table S4).
To place the samples analyzed in a phylogenetic context, 215 GenBank sequences obtained from species distributed among different groups within Rickettsia were downloaded and used for phylogenetic analyses (Table S5).Note that, aiming to cover the entire diversity of the genus, a greater number of GenBank sequences were initially included in exploratory analyses.However, long-branch attraction artefacts were identified (for example, for the sequences of R. bellii and R. felis), so they were not included in the final datasets.The sequences ultimately selected corresponded to 120 isolates for which ompA and gltA gene data were available, including 40 isolates designated "type strain", which were used to delineate species within the genus Rickettsia, and 16 undetermined or uncultured Rickettsia species.
Homologous GenBank and newly obtained sequences were automatically aligned using MAFFT with the E-INS-i algorithm and default settings, as implemented in Geneious.Low-quality ends were trimmed before the alignment errors were manually adjusted.Both alignments used in this study are available as Supporting Information (Alignments S1 and S2).
For each gene, two common phylogenetic approaches, i.e., maximum likelihood (ML) and Bayesian inference (BI), were used to determine the phylogenetic position of the isolates obtained within the genus Rickettsia.We used different phylogenetic tools available in CIPRES (CyberInfrastructure for Phylogenetic Research; www.phylo.org).Specifically, single-gene ML trees were inferred using IQ-TREE v. 2.1.2[22].The bestfit model of nucleotide substitution and the optimal partitioning scheme for each gene (both were divided into three partitions corresponding to each codon position) were selected by the integrated version of ModelFinder [23].The complete bootstrap option with 1000 nonparametric bootstrap replicates was used to determine bootstrap support (BS) for each node.
The BI individual analyses were performed using the Metropolis-coupled Markov chain Monte Carlo (MCMCMC) method, as implemented in MrBayes v. 3.2.7 [24].For each gene dataset, the best-fit substitution model was estimated by sampling across the substitution model space with the reversible-jump Markov chain Monte Carlo (MCMC) method [25], allowing gamma distributed rate heterogeneity across sites and a proportion of invariant sites.The partitioning scheme selected by ModelFinder was used, unlinking model parameters across different partitions.Four independent runs of 20 million generations, each with six chains, were used, sampling trees every 1000 generations, with the first 25% being discarded as burn-in.The posterior probabilities (PPs) were calculated from the remaining probabilities.Bayesian analyses were automatically stopped when the average standard deviation of the split frequencies was smaller than 0.01.In addition, the convergence of the runs was checked in Tracer v. 1.7.2 [26] by checking that the effective sampling size (ESS) values for all parameters fell below 200.Single-gene trees were visually inspected in FigTree v. 1.4.4 [27] before being processed with Adobe Illustrator CS5 (Adobe Systems Inc., San Jose, CA, USA).

Statistical and Graphical Analyses
Fisher's one-way repeated measures ANOVA (parametric) was used to estimate the difference in means of infection levels among tick species, life stages, months, and years.A post hoc test was performed with a Holm p-value correction due to multiple comparisons.A significant difference was assumed for values of p < 0.05.We estimated the effect size using Hedges' G with a 95% confidence interval (CI).For the Bayesian interpretation, a Cauchy distribution r = 0.707 was considered a priory, the highest density interval (HDI) of 95% as credible intervals, and an interpretation of the logarithm of the Bayes factor according to the interpretation scales suggested by Jeffreys [28].All the analyses were performed using the package ggstatsplot [29] in R [30].
To obtain a visual representation of the relationship between ticks and their associated Rickettsia, a heat map was made.First, we loaded the frequency of infected ticks of different species and converted these data into a matrix using R, and then the pheatmap function of the pheatmap package [31] was used to build the map.

Tick Identification
A total of 7397 ticks removed from 7180 people were identified at the species level, with only Rh. sanguineus (s.l.) ticks being identified at the group level.The tick specimens identified included 140 larvae, 2227 nymphs, and 5030 adults belonging to five genera and eight species (Table 1).Ixodes ricinus was the most represented (427/3269, 44.19%).In general, the highest number of tick bites was reported during spring and summer (April to July), with the peak occurring in May and June of most years.However, from 2021 onwards, the activity peak began earlier, in April, and remained very high until July.

Molecular Detection of Rickettsia
Rickettsial DNA was detected in 1177 of the 7397 ticks analyzed (Table 1).For these specimens, gltA, ompA, or both genes could be amplified, which represents an overall prevalence of 15.91%.Although I. ricinus was the species most frequently carrying Rickettsia, the highest percentage of infected specimens per species corresponded to the genera Dermacentor (D. marginatus, 35.88% and D. reticulatus, 20.19%), followed by Rhipicephalus (Rh.sanguineus, 19.5% and Rh.bursa, 16%) (Table 1).In terms of the percentage of ticks infected according to their life stage and sex, we found that 80% of the positive ticks were adults (54% females and 26% males), 1% were larvae, and 19% were nymphs.In terms of monthly differences, ticks infected by different Rickettsia species were detected throughout the year (Figure 2), with June being the month with the highest percentage of Rickettsia-positive ticks (24.06%), followed by May (19.73%) and April (12.8%).The interannual distribution of Rickettsia-carrying ticks was quite regular throughout the study period, showing a decreasing trend.the year (Figure 2), with June being the month with the highest percentage of Rickettsiapositive ticks (24.06%), followed by May (19.73%) and April (12.8%).The interannual distribution of Rickettsia-carrying ticks was quite regular throughout the study period, showing a decreasing trend.

Sequences and BLAST (Identification of Rickettsia Species)
A total of 504 isolates out of the 1177 Rickettsia-positive samples were sequenced, resulting in 581 partial sequences (311 sequences corresponding to ompA and 270 to gltA, Table S4), which allowed the identification of 10 different taxa.Specifically, based on BLAST searches of the ompA gene sequences, all our tick-isolated Rickettsia sequences showed high similarity to GenBank sequences corresponding to nine different members of the spotted fever group I (SFGI): R. aeschlimannii, R. conorii subsp.conorii, R. conorii subsp.raoultii, R. massiliae, R. monacensis, R. sibirica subsp.mongolitimonae, R. slovaca, Candidatus R. barbariae, and Ca.R. rioja.Similarly, most isolates positive for the gltA gene matched with some Rickettsia species belonging to the SGFI, and 15 isolates matched with the species R. helvetica.

Sequences and BLAST (Identification of Rickettsia Species)
A total of 504 isolates out of the 1177 Rickettsia-positive samples were sequenced, resulting in 581 partial sequences (311 sequences corresponding to ompA and 270 to gltA, Table S4), which allowed the identification of 10 different taxa.Specifically, based on BLAST searches of the ompA gene sequences, all our tick-isolated Rickettsia sequences showed high similarity to GenBank sequences corresponding to nine different members of the spotted fever group I (SFGI): R. aeschlimannii, R. conorii subsp.conorii, R. conorii subsp.raoultii, R. massiliae, R. monacensis, R. sibirica subsp.mongolitimonae, R. slovaca, Candidatus R. barbariae, and Ca.R. rioja.Similarly, most isolates positive for the gltA gene matched with some Rickettsia species belonging to the SGFI, and 15 isolates matched with the species R. helvetica.
According to BLAST analyses (Table S3), the most prevalent Rickettsia species was R. massiliae (24%, 121/504).It was detected mainly in ticks of the genus Rhipicephalus (Rh.bursa, 67/121 and Rh.sanguineus, 45/121) (Figure 3).The second most prevalent species was R. monacensis, with a prevalence of 22.22% (Table S3, 112/504), which was mainly detected in I. ricinus (106/111).These were followed by a group consisting of R. slovaca, R. aeschlimanii, and R. conorii subsp.raoultii, with prevalences of 16.43%, 15.63%, and 13.82%, respectively.Rickettsia slovaca and R. conorii subsp.raoultii were mainly detected in D. marginatus, and Hy.marginatum and R. aeschmimanii were detected in Hy. marginatum.In the case of Ca.R. rioja and R. helvetica, the prevalence rates were very similar (3.6% and 2.8%, respectively), showing in both cases a high specificity for the vector species: Ca.R. rioja was almost exclusively found in D. marginatus specimens, and R. helvetica was found in I. ricinus.For the remaining members of Rickettsia found in this study (R. conorii subsp.conorii, Ca.R. barbariae, and R. sibirica subsp.mongolitimonae), their prevalence ranged from 0.2 to 0.4%, with each being found in only one tick species: R. conorii subsp.conorii and Ca.R. barbariae in Rh. bursa, and R. sibirica subsp.mongolitimonae in Hy. marginatum.

Phylogenetic Analyses
After excluding primer binding regions and low-quality ends, the alignments used for phylogenetic purposes were 608 and 341 bp long (ompA and gltA, respectively).The optimal partitioning scheme and best-fit model of nucleotide substitution for each gene matrix used in our phylogenetic analyses are shown in Table S6.For both genes, both ML and BI analyses yielded similar topologies, so only the Bayesian tree showing PP and BS values is shown (Figures 4 and 5).For better visualization, well-supported clades comprising a high number of isolates were collapsed (solid black triangles).The extended versions of these trees can be found in the supplementary material (Figures S1 and S2).For the remaining members of Rickettsia found in this study (R. conorii subsp.conorii, Ca.R. barbariae, and R. sibirica subsp.mongolitimonae), their prevalence ranged from 0.2 to 0.4%, with each being found in only one tick species: R. conorii subsp.conorii and Ca.R. barbariae in Rh. bursa, and R. sibirica subsp.mongolitimonae in Hy. marginatum.

Phylogenetic Analyses
After excluding primer binding regions and low-quality ends, the alignments used for phylogenetic purposes were 608 and 341 bp long (ompA and gltA, respectively).The optimal partitioning scheme and best-fit model of nucleotide substitution for each gene matrix used in our phylogenetic analyses are shown in Table S6.For both genes, both ML and BI analyses yielded similar topologies, so only the Bayesian tree showing PP and BS values is shown (Figures 4 and 5).For better visualization, well-supported clades comprising a high number of isolates were collapsed (solid black triangles).The extended versions of these trees can be found in the Supplementary Material (Figures S1 and S2).
Individually, none of the two molecular regions analyzed confidently resolved all relationships among the newly obtained isolates and different Rickettsia species used as references, as indicated by the relatively low support values recovered for some internal branches in both single-gene trees (Figures 4, 5, S1 and S2).This is especially true for the gltA tree (Figures 5 and S2), with many isolates arising from a polytomy and containing sequences obtained from the type strains of several other species.However, in general, most of our isolates formed different well-defined species clusters in both the ompA and gltA trees.An interrupted branch (//) indicates its length has been reduced.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.An interrupted branch (//) indicates its length has been reduced.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.An interrupted branch (//) indicates its length has been reduced.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.An interrupted branch (//) indicates its length has been reduced.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.
The phylogenetic tree inferred from our gltA partial sequences (Figures 5 and S2) is highly unresolved, as evidenced by the existence of a large poorly supported polytomy (PP = 0.74, BS = 76%), corresponding to the SFGI.The other main groups within the genus Rickettsia are either well or moderately supported, except for SFGII, which was recovered as polyphyletic.Specifically, the CG and the group formed by R. helvetica and R. asiatica were monophyletic with high support (PP = 1, BS = 95% and PP = 1, BS = 89%, respectively), the TIG was only moderately supported (PP = 1, BS = 63%), the typhus group (TG) was fully supported (PP = 1, BS = 100%), and the sister had low support for the polytomy (PP = 1, BS = 47%).Most of our isolates (207 out of 270) form part of the mentioned polytomy, with only some being included in the subgroups arising from it.In detail, one isolate (Isolate 420) clustered within the "Ca.R. barbariae" group, only moderately supported by ML (PP = 0.85, BS = 64%); three isolates identified as Rickettsia conorii subsp.raoultii by BLAST and two undetermined ones grouped together in a well-supported cluster (PP = 0.97, BS = 78%) that does not comprise any sequences from a type strain and, thus, is unnamed; 51 isolates form an unsupported group (PP = 0.55, BS = 15%) that includes the type strain of R. slovaca; and 72 additional ones clustered together with the type strain of R. massiliae, also forming an unsupported group (PP = 0.75, BS = 57%).The remaining isolates were distributed in different subgroups not arising from the polytomy as follows: 15 isolates clustered within the highly supported clade "R.helvetica + R. asiatica" (PP = 1, BS = 89%), and 48 isolates grouped together with the type strain of R. monacensis (PP = 0.95, BS = 54%).

Distribution of Rickettsia Species within the Sampling Area
The geographical distribution of some Rickettsia species is clearly delimited within Castilla y León (Figure 6).This is the case for Ca.R. rioja, since 100% of the ticks infected by this bacterium were found in the northwestern region.Additionally, in the northern region, we found 66.7% of all ticks that carry R. slovaca in Castilla y León.
Most R. aeschlimanii (60%), R. helvetica (71.39%),R. massiliae (47.9%), and R. monacensis (65.5%) isolates were recovered from ticks found in the southern part of the sampling area.In the case of the ticks infected by R. conorii subsp.raoultii, however, they were more evenly distributed throughout the study area.In turn, most ticks infected by R. conorii subsp.conorii were found in the north, and those infected by Ca.R. barbariae and R. sibirica subsp.mongolitimonae were found in the south, although these infections were very infrequent.
(65.5%) isolates were recovered from ticks found in the southern part of the sampling area.In the case of the ticks infected by R. conorii subsp.raoultii, however, they were more evenly distributed throughout the study area.In turn, most ticks infected by R. conorii subsp.conorii were found in the north, and those infected by Ca.R. barbariae and R. sibirica subsp.mongolitimonae were found in the south, although these infections were very infrequent.

Statistical Analyses
Fisher's one-way repeated measures ANOVA test revealed that, considering only the infected ticks, there were not statistically significant differences among different years (Figure 7A).study).(B) Sampling between 2018 and 2022 (this study).

Statistical Analyses
Fisher's one-way repeated measures ANOVA test revealed that, considering only the infected ticks, there were not statistically significant differences among different years (Figure 7A).Although we originally observed statistically significant differences among months, after post hoc pairwise t-tests with the Holm p-value correction, no significant differences existed (Figure 7B).
There were no statistically significant differences among different stages (Figure 8A), although they were observed among tick species (Figure 8B).The effect size (ωp = 0.82) was large, as per Field's conventions [32].The Bayes Factor for the same analysis revealed that the data were more probably under the alternative hypothesis, as compared to the (B) (A) Although we originally observed statistically significant differences among months, after post hoc pairwise t-tests with the Holm p-value correction, no significant differences existed (Figure 7B).
There were no statistically significant differences among different stages (Figure 8A), although they were observed among tick species (Figure 8B).The effect size (ωp = 0.82) was large, as per Field's conventions [32].The Bayes Factor for the same analysis revealed that the data were more probably under the alternative hypothesis, as compared to the null hypothesis (BF01-22.81).This can be considered very strong evidence in favor of the alternative hypothesis [28].This global effect was carried out by post hoc pairwise t-tests, which revealed that there are significant differences among different species (Figure 8).
null hypothesis (BF01-22.81).This can be considered very strong evidence in favor of the alternative hypothesis [28].This global effect was carried out by post hoc pairwise t-tests, which revealed that there are significant differences among different species (Figure 8).

Discussion
Ticks and tick-borne diseases (TBDs), especially tick-borne rickettsioses (TBR), are a growing problem for human health and have been reported in different countries and regions, including neighboring countries to Spain, i.e., France and Portugal [33,34], among many others.Moreover, the incidence of emerging TBD, including rickettsiosis, is likely to increase in the near future [35,36].For these reasons, studies addressing the surveillance of ticks and emerging and re-emerging TBDs, such as rickettsial diseases, are essential not only for identifying areas and populations at risk but also for implementing appropriate prevention and control measures.These prevention measures include using appropriate clothes and repellent against ticks when visiting risk areas, self-examination

Discussion
Ticks and tick-borne diseases (TBDs), especially tick-borne rickettsioses (TBR), are a growing problem for human health and have been reported in different countries and regions, including neighboring countries to Spain, i.e., France and Portugal [33,34], among many others.Moreover, the incidence of emerging TBD, including rickettsiosis, is likely to increase in the near future [35,36].For these reasons, studies addressing the surveillance of ticks and emerging and re-emerging TBDs, such as rickettsial diseases, are essential not only for identifying areas and populations at risk but also for implementing appropriate prevention and control measures.These prevention measures include using appropriate clothes and repellent against ticks when visiting risk areas, self-examination after outdoor activities, the removal of any ticks as soon as possible, and the deparasitation of domestic and farm animals.
The present study provides relevant information on the diversity and seasonal distribution of human-biting ticks in Castilla y León (northwest Spain) over a five-year span, as well as relevant molecular data on the tick-associated Rickettsia that they could transmit to humans in this region.Compared to what was observed by Fernández-Soto [14], the most evident difference is the changing distribution pattern of two species, i.e., R. slovaca and R. monacensis, within the study area.We have seen here that R. slovaca is now restricted to the north of Castilla y León, especially to the northwest region, while years ago, it was found in virtually the whole study area.However, in the case of R. monacensis, we have seen the opposite trend: previously, this species was almost restricted to a single locality in the south, while it is now found across Castilla y León, although it is still more prevalent in southern areas.It is also important to note that R. monacensis is the species that increased the most in number and distribution area.
We analyzed both the ompA and gltA benchmark genes and found that our isolates belong to different members of the SFGI and to the species R. helvetica, which is currently unassigned to a defined group within the genus.
Specifically, based on the information provided by ompA (Figures 4 and S1), most of our isolates clustered within the SFGI and corresponded to Ca. R. rioja, Ca.R. barbarie, R. aeschlimannii, R. conorii subsp.raoultii, and R. massiliae.R. monacensis, and R. slovaca.Considering the phylogeny inferred from the gltA gene (Figures 5 and S2), which is not resolved as evidenced by the polytomy of indeterminate relatedness, the existence of different independent well-supported clades adds some information to that provided by ompA: all 15 isolates identified as R. helvetica by BLAST form a group sister to R. asiatica with high support, as previously found using phylogenomic data [37].Moreover, 48 gltA isolates belong to R. monacensis.
Focusing on the ticks biting people in the study area, Ixodes ricinus was clearly the predominant species (it was the one most frequently removed from people), followed by Rh. bursa, which coincides with the results obtained in other studies carried out by our group [13].However, we detected an increase in the percentage of three tick species: Hy. marginatum (11.81% vs. 9.08%), Hy. lusitanicum (11.81 vs. 9.08%), and Rh.bursa (17.56% vs. 12.15%).This upwards trend was already noted during the last years of our previous study and has been maintained over time.For the remaining tick species, there were no significant changes in terms of abundance.
Notably, 15% of the ticks tested were positive for rickettsial DNA, a much greater percentage than that observed in a previous study (approximately 7%) carried out by our group in the same area years ago [14], but very similar to those more recently reported in some locations of Italy (15-17%) [38,39] and Poland [40].This percentage is, however, much lower than those reported in other European countries, such as France [34], the Netherlands [41], Serbia [42], and Turkey [43], ranging from 22 to 44%.In contrast, in Sweden [44] and Romania [45], the prevalence rates are less than 9%.
We found that Dermacentor marginatus was the species with the highest prevalence of rickettsial infection (35.88%), which is consistent with what was observed in the same area by Fernández-Soto [14] but also in other regions of Spain [46,47], Italy [38], and France [34].
Seven members of the genus Rickettsia, i.e., Rickettsia aeschlimannii, R. conorii subsp.conorii, R. conorii subsp.raoultii, R. massiliae, R. monacensis, R. helvetica and R. slovaca, have been found to be associated with different species of ticks in this study but also in that previously conducted by Fernández-Soto [14] in the same area, which clearly indicates that all these bacteria were present in Castilla y León since at least 20 years ago.Nevertheless, we also found several Rickettsia species that were nonexistent, or at least not identified, in the mentioned study.This is the case for Candidatus R. rioja, Candidatus R. barbariae, and R. sibirica subsp.mongolitimonae.However, it is important to note that these species have also been found in regions close to the sampling area.Specifically, Candidatus R. rioja, an SFG Rickettsia species involved in DEBONEL/TIBOLA cases, has been detected in patients bitten by ticks in La Rioja, a region to the northeast of Castilla y León [48].Similarly, R. sibirica subsp.mongolitimonae has been found in Rh. bursa removed from cattle and Rh.pusillus from rabbits in central Spain [49].Most importantly, this species has been reported to cause human infections [33,50].In regions more distant from the sampling area, Candidatus R. barbariae, originally named strain PoTiRb169, has also been detected in Rh. bursa ticks collected in Sardinia [51].These findings, along with ours, confirm the broad distribution of these agents in southern Europe but also raise the possibility that some cases of rickettsiosis diagnosed as MSF, mostly caused by R. conorii subsp.conorii, could actually be due to other members of Rickettsia.
By comparing the Rickettsia species currently present in the study area with those previously reported, it can be seen that some of them have increased in number to the detriment of others that have decreased.In both periods (1997-2002 and 2018-2022), the most prevalent Rickettsia species in Castilla y León was R. massiliae, but curiously, there were differences in the tick vector species.In the study carried out between 1997 and 2002, all the ticks carrying R. massiliae belonged to the Rh.sanguineus complex, while, in our study, they belonged not only to the Rh.sanguineus complex but also to the species Rh. bursa.We detected more ticks infected by this Rickettsia species in the southern region of the community, which was expected considering that its main vectors are mainly found in this area; however, compared to what was reported in a previous study, the greatest increase in Rickettsia occurred in the northern provinces.This could be explained by the fact that both Rh.bursa and Rh.sanguineus are beginning to establish themselves in the northern provinces of the community [13].
According to our results, currently, the second most prevalent Rickettsia species is R. monacensis, which is somewhat different from what was reported in the study by Fernández-Soto [14], in which R. slovaca and R. aeschlimanii were more abundant.Rickettsia monacensis, mainly detected in I. ricinus, is the species that increased the most in terms of number and distribution range with respect to the mentioned study.Although it is still more abundant in the south, we have seen an increase in the north (northeast), which coincides with the increase in I. ricinus in this area.
The prevalence of ticks infected by R. slovaca and R. aeschlimanii has decreased compared with that found in a previous study, while the number of ticks infected by R. conorii subsp.raoultii has increased over time (Table S7).The explanation for this increase is that R. conorii subsp.raoultii shares vectors with the other two mentioned Rickettsia species.According to previous results [14], R. conorii subsp.raoultii was mostly transmitted by ticks of the genus Dermacentor, while our data indicate that it is now transmitted by ticks belonging to Dermacentor, but also to the genus Hyalomma in similar proportions, as depicted in Figure 3.
It is very important to highlight the role of Hy. marginatum as a vector of this Rickettsia species since, in a previous study, it was mainly transmitted by the genus Dermacentor.For the other two Rickettsia species, ticks of the genera Dermacentor (mainly D. marginatus) and Hyalomma (Hy.marginatum) were the main vectors of R. slovaca and R. aeschlimanii, respectively.
When studying the distribution of these Rickettsia species, we observed that R. slovaca continues to be predominant in the north, although it has moved to the west, as well as its main vector, D. marginatus [13], and that R. aeschlimanii has a very similar distribution to that seen years ago, with a higher prevalence in the south and an absence in the center of the community.Rickettsia helvetica was less prevalent in the current study than in that by Fernández-Soto [14], although in both studies, it was found almost exclusively in I. ricinus and mostly in the southern region.
Ca. R. rioja was detected for the first time in ticks attached to humans in Castilla y León in 2019.It should be noted that this uncultivated species, which is phylogenetically very close to R. conorii subsp.raoultii, can only be distinguished from the latter by analysis of ompA because most authors include it under R. conorii subsp.raoultii [5].However, since 2019, the incidence of Candidatus R. rioja has increased until 2022, when the highest number of specimens of this genus of Rickettsia was detected.This species is confined to the northern zone, although, as we have seen for R. slovaca, it is moving westwards.This geographic shift seems logical since both Candidatus R. rioja and R. slovaca share the same vector.
The species with the greatest number of ticks infected by Rickettsia species was Rh. bursa (eight species), followed by I. ricinus (six), D. marginatus, and Hy.marginatum (five).These results are quite different from those observed by Fernández-Soto [14], who identified I. ricinus as the main vector species associated with six different Rickettsia species.
The results of this study show important aspects related to human-biting ticks and their associated pathogenic Rickettsia.However, this is only a partial vision of the human-animalenvironment interface.Considering that the health of humans, animals, and ecosystems are interconnected, it is urgently needed to obtain a whole picture of the actual situation, including the most medically/veterinary important pathogens transmitted by ticks.Future directions include multidisciplinary collaborations and cross-sectoral studies for a better understanding of this multifaceted problem and to be able to tackle the potential emergence of tick-related zoonotic diseases.

Conclusions
This long-term study represents a step forward to a better understanding of the Rickettsia associated with ticks in northwestern Spain.Despite the fact that it was impossible to molecularly characterized all ticks positive for Rickettsia, our comprehensive morphological study of the ticks and the identification of a considerable number of bacterial isolates, at molecular level, allow us to have a clear vision of the epidemiological situation in the region.
Our results show the presence of a variety of pathogenic Rickettsia in ticks removed from humans in Castilla y León, including several species that have gone undetected to date.We have observed a wider geographical and seasonal expansion of tick populations than previously acknowledged in this area, as well as an increase in the prevalence of ticks carrying Rickettsia, with a consequent risk to human health.
When comparing our results with those obtained in previous studies conducted in the same area, evident changes in the distribution patterns of different rickettsial species were observed.Ixodes ricinus and Rh.bursa are the predominant tick species removed from humans.Ten different members of the genus Rickettsia (nine SFGI and R. helvetica) were detected in the ticks collected in this study, with R. massiliae being the most prevalent species and R. monacensis being the species with the greatest increase in number and distribution area.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects15080571/s1, Figure S1: Bayesian phylogenetic tree of the rickettsial isolates obtained, based on 608 nucleotide positions of ompA; Figure S2: Bayesian phylogenetic tree of the rickettsial isolates obtained, based on 341 nucleotide positions of gltA; Table S1: Primers used for PCR amplification and size of the amplicons; Table S2: Cycling parameters used in this study; Table S3: Rickettsial isolates obtained from different tick species, and putative identity based on BLAST results of ompA and gltA partial sequences; Table S4: Rickettsial isolates newly obtained from different tick species in this study; Table S5: GenBank sequences generated from 120 strains, included in the phylogenetic analyses conducted in this study; Table S6: Best partition scheme and substitution model(s) for both datasets used in phylogenetic analyses, chosen according to BIC; Table S7: Comparison of the prevalence of three Rickettsia species in different periods in Castilla y León.Alignment S1: alignment of ompA gene sequences; Alignment S2: alignment of gltA gene sequences.References in Supplementary Materials: [7,14,20,21,37].

Figure 2 .
Figure 2. Number of Rickettsia-positive and Rickettsia-negative ticks per species versus months.

Figure 2 .
Figure 2. Number of Rickettsia-positive and Rickettsia-negative ticks per species versus months.

Figure 3 .
Figure 3. Heat map showing the number of Rickettsia isolates per tick species.Red represents the highest number of infected ticks and yellowish the smallest number.

Figure 3 .
Figure 3. Heat map showing the number of Rickettsia isolates per tick species.Red represents the highest number of infected ticks and yellowish the smallest number.

Figure 4 .
Figure 4. Bayesian phylogenetic tree of the rickettsial isolates obtained, based on 608 nucleotide positions of ompA, with members of the "Canadensis group" as outgroup.Bayesian posterior probabilities (PP) and maximum likelihood bootstrap support (BS) values are shown above and below each branch if PP > 0.90 and BS > 50%, respectively.Lower support values are represented by dashes.Solid dots indicate nodes fully supported (PP = 1, BS = 100%).For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.An interrupted branch (//) indicates its length has been reduced.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.

Figure 4 .
Figure 4. Bayesian phylogenetic tree of the rickettsial isolates obtained, based on 608 nucleotide positions of ompA, with members of the "Canadensis group" as outgroup.Bayesian posterior probabilities (PP) and maximum likelihood bootstrap support (BS) values are shown above and below each branch if PP > 0.90 and BS > 50%, respectively.Lower support values are represented by dashes.Solid dots indicate nodes fully supported (PP = 1, BS = 100%).For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.An interrupted branch (//) indicates its length has been reduced.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.

Figure 5 .
Figure 5. Bayesian phylogenetic tree of the rickettsial isolates obtained, based on 341 nucleotide positions of gltA, with members of the "Canadensis group" as outgroup.Bayesian posterior probabilities (PP) and maximum likelihood bootstrap support (BS) values are shown above and below each branch if PP > 0.90 and BS > 50%, respectively.Lower support values are represented by dashes.Solid dots indicate nodes fully supported (PP = 1, BS = 100%).For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.An interrupted branch (//) indicates its length has been reduced.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.

Figure 5 .
Figure 5. Bayesian phylogenetic tree of the rickettsial isolates obtained, based on 341 nucleotide positions of gltA, with members of the "Canadensis group" as outgroup.Bayesian posterior probabilities (PP) and maximum likelihood bootstrap support (BS) values are shown above and below each branch if PP > 0.90 and BS > 50%, respectively.Lower support values are represented by dashes.Solid dots indicate nodes fully supported (PP = 1, BS = 100%).For newly obtained sequences, the isolate number is followed by the species name.In the case of the GenBank sequences used as reference, the species name is followed by strain ID and accession number, with type strains appearing in bold and marked with a (T).For representation purposes, well-supported clades comprising isolates obtained here plus one type strain have been collapsed, thus appearing as solid black triangles.Continuous and discontinuous vertical lines represent mono-and paraphyletic groups, respectively.The scale bar represents the average number of substitutions per site.

Figure 6 .Figure 6 .
Figure 6.Distribution of Rickettsia species within the sampling area (Castilla y León, northwestern Spain), in two different periods.Each species is represented by a colored circle with size Figure 6.Distribution of Rickettsia species within the sampling area (Castilla y León, northwestern Spain), in two different periods.Each species is represented by a colored circle with size proportional to the number of isolates found.(A) Sampling period between 1997 and 2002 (previous study).(B) Sampling between 2018 and 2022 (this study).

Figure 7 .
Figure 7. Differences in the number of infected ticks per year (A) and month (B).

Figure 7 .
Figure 7. Differences in the number of infected ticks per year (A) and month (B).

Figure 8 .
Figure 8. Differences in the number of infected ticks per stages (A) and species (B).

Figure 8 .
Figure 8. Differences in the number of infected ticks per stages (A) and species (B).

Table 1 .
Data on the ticks analyzed, including tick species, number of specimens collected for each species, and number of specimens carrying Rickettsia species.